rm(list =ls())
options(scipen=999)
gc()
packages <-c("tidyverse","estimatr","plm","stargazer",
             "fastDummies","ICCbin","ihs","readstata13","xtable",
             "arm")

new.packages <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

lapply(packages, require, character.only = TRUE)
rm(packages, new.packages)

setwd("PUT YOUR DIRECTORY HERE")

## insheet data
data <- read.csv("./Datasets/panel_bare_bones.csv")

## clean data
data$X <- NULL
data$locality_year <- NULL
data$subbotnik_inkind <- NULL
data$locality <- as.character(tolower(data$locality))

# insheet population 
dataa <- read.csv("./Datasets/main_confounders.csv")
data <- merge(data, dataa[,c("pop", "locality")], by = c("locality"), all.x = T)
rm(dataa)


###############################
### REGRESSIONS FOR STRIKES ###
###############################

#subset to years with observations:
data_t <- data[which(data$year>1979),]

#### OLS (averaging over years)
data_agg <- aggregate(data_t[,c("strikes", "commanders", "pop")], by=list(Category=data_t$locality_numeric), FUN=mean, na.rm=T)

## create normalized by population
data_agg$commanders_pop <- data_agg$commanders/data_agg$pop

## run model
m_strikes_ols_averaged_pop <- lm(scale(strikes) ~ scale(commanders_pop), data = data_agg)


################################
### REGRESSIONS FOR SABOTAGE ###
################################

# subset to years with observations:
data_t <- data[which(data$year<1980),]

# complete cases
data_t <- data_t[,c("locality_numeric", "commanders", "subbotnik_inkind_", "pop")]
data_t <- na.omit(data_t) 

# turn subbotnik upside down to mean sabotage
data_t$subbotnik_inkind_ <- data_t$subbotnik_inkind_*-1

#### OLS (averaging over years)
data_agg <- aggregate(data_t[,c("subbotnik_inkind_", "commanders", "pop")], by=list(Category=data_t$locality_numeric), FUN=mean, na.rm=T)

## create normalized by population
data_agg$commanders_pop <- data_agg$commanders/data_agg$pop

## put out models
m_subbotnik_inkind__ols_averaged_pop <- lm(scale(subbotnik_inkind_) ~ scale(commanders_pop), data = data_agg)



## save as tex
stargazer(m_strikes_ols_averaged_pop, m_subbotnik_inkind__ols_averaged_pop, 
          dep.var.labels = c("Protest (z)", "Sabotage (z)"),
          style = "qje",
          covariate.labels = c("Repression (z)"),
          star.char = c("*", "**", "***"),
          star.cutoffs = c(0.1, 0.05, 0.01),
          omit = c("Constant"),
          omit.stat = c("rsq", "f", "ser"), 
          omit.table.layout = "n",
          add.lines = list(c("Fixed effects", "No", "No"),
                           c("Controls", "No", "No")),
          out = "PUT YOUR FILEPATH HERE")

